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Abstract. A number of supernova remnants (SNRs) show nonthermal X-rays assumed to be synchrotron emission 
from shock accelerated TeV electrons. The existence of these TeV electrons strongly suggests that the shocks in 
SNRs are sources of galactic cosmic rays (CRs). In addition, there is convincing evidence from broad-band studies 
of individual SNRs and elsewhere that the particle acceleration process in SNRs can be efficient and nonlinear. If 
SNR shocks are efficient particle accelerators, the production of CRs impacts the thermal properties of the shock 
heated, X-ray emitting gas and the SNR evolution. We report on a technique that couples nonlinear diffusive 
shock acceleration, including the backreaction of the accelerated particles on the structure of the forward and 
reverse shocks, with a hydrodynamic simulation of SNR evolution. Compared to models which ignore CRs, the 
most important hydrodynamical effects of placing a significant fraction of shock energy into CRs are larger shock 
compression ratios and lower temperatures in the shocked gas. We compare our results, which use an approximate 
description of the acceleration pro cess, with a more complete model where the full CR transport equations are 
solved (i.e., iBerezhko et al.L l2002il . and find excellent agreement for the CR spectrum summed over the SNR 
lifetime and the evolving shock compression ratio. The importance of the coupling between particle acceleration 
and SNR dynamics for the interpretation of broad-band continuum and thermal X-ray observations is discussed. 
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1. Introduction 

It is commonly believed that the shocks in super- 
nova remnants (SNRs) produce the majority of galac- 
tic cosmic rays (CRs) with energies b elow ~ 10 15 
eV via diffusive shock ac celeration (e.g., IPrurvL Il983_ 
iBlandford fc Eichlerl Il987l) . Convincing support for the 
production of TeV electrons in SNRs comes from the 
sync hrotron interpret ation of nonthermal X-ray emission 
(e.g.. iRevnoldsl Il998j) observ ed in an increas i ng nu mber of 
young SNRs such a s SN1006 llKovama et allll995lh C as A 
.Allen et al.lll997lh G347.3-0.5 llSlane et allllQQfllh and 
RCW 86 llBorkowski et aflkOOlh . 1 
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1 As of this writing, there is no una mbiguous evidence for 
the production of TeV ions in SNRs fsee lBerezhko et all 120021 
for a discussion of SN1006 in this regard). The recent claim 
that TeV emission from SNR RX J1713. 7-3946 (also called 
G347 3-0.5) observed by C ANGAROO II is from pion-decay 
(e.g., |Ejiarnoto_^t_alJ, |2002j) i s still under deba te (see, for ex- 
ample. iReimer fc Pohj 12002. iButt et all I2002T) . 



In addition to their putative role in accelerating cos- 
mic rays, the shocks in SNRs heat the ambient interstel- 
lar medium and ejecta to X-ray emitting temperatures. 
The interpretation of these X-ray observations leads to 
inferences for important quantities such as the supernova 
(SN) explosion energy, ejecta mass and composition, am- 
bient densities, shock speed, and rate of electron and 
proton equilibration. That CR production and thermal 
heating in SNRs may be coupled comes from the fact 
that diffusive shock acceleration is intrinsically efficient in 
high Mach number shocks if even a small fraction of the 
shock hea ted plasma is injected into the acceleration pro- 
cess (e.g-lEllison fc Eichlerl ll984Ujones fc Ellisorl Il991_ 
IBerezhko fc Ellisorlll999UMalkovMl998HBlasil2002|h 

If strong coupling between particle acceleration and 
shock heating occurs, the modeling of efficient particle ac- 
celeration in SNRs offers the possibility of using the high- 
quality X-ray and 7-ray data currently being collected by 
spacecraft (e.g., Chandra, XMM-Newton, INTEGRAL) to 
address fundamental questions concerning the SNR origin 
of CRs and the underlying physics of diffusive shock accel- 
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eration, particularly the injection of thermal particles into 
the acceleration process. Furthermore, strong coupling im- 
plies that the inferences made from X-ray observations 
may differ substantially between interpretations which in- 
clude particle acceleration self-consistently and those that 
do not. 

Despite the expected efficiency of diffusive shock 
acceleration, X-ray spectra from SNRs have generally 
been modeled and interpreted assuming that the shocks 
place an insignificant fraction of their energy in cos- 
mic rays. Excep ti ons to th i s inclu de th e early works of 
IChevalierl lll983l) iHeayend l|l984|) . and iBoulares fc Col 
<ll988y Chevalier! Jl983l) investigated the effects of cosmic- 
ray pressure on SNR dynamics using a two-fluid, self- 
similar solution with an arbitrary fraction of thermal gas 
(adiabatic index 7 = 5/3) and relativi stic gas (7 = 4/3). 
More recen t work has b een done by iDorfi fc Bohringerl 
ill 9931) and iDorfil ll 19941) . In our preliminary work (i.e., 
iDecourchelle et all I2OO0I) . we developed a model which 
coupled the approximate nonlin e ar (N L) acceleration cal- 
culation of iBerezhko fe Ellisonl lll999l) with an analytic, 
self-similar desc ription of the SNR hydrodynamics (i.e., 
IChevalierl Fl983T) and a non- equilibrium ionization calcula - 
tion of X-ray emission ( e.g.. IDecourchelle fc Balletl ll 994h . 
We illustrated this model by fitting ASCA and RXTE ob- 
servations of Kepler's SNR and found adequate fits when 
acceleration was efficient at the forward shock but ineffi- 
cient at the reverse shock. 

The effects of efficient particle acceleration on SNR hy- 
drodynamics where calculated in a hydro com puter sim- 
ulation of SNRs bv iBlondin fc Ellisonl (|200ll) . This was 
done by globally changing the effective ratio of specific 
heats, 7 c ff, from 5/3 to values approaching 1, but did not 
include coupling between the acceleration and the hydro. 
As 7off was decreased, the shocked gas became more com- 
pressible, the shock compression ratio increased, and the 
interactio n region between the forwa rd and reverse shocks 
narrowed. IBlondin fc Ellisonl were able to show in 

two and three-dimensional simulations that if the inter- 
action region was narrow enough, convective instabilities 
produced Rayleigh- Taylor fingers of dense ejecta material 
which were able to reach and perturb the forward shock. 

Here, we introduce and describe in detail a CR-Hydro 
mod el which uses the same NL acceleration calculation 
of lBerezhko fc Ellisonl lll999|). but replaces th e self-similar 
description used in IDecourchelle et alJ ll2000l) with a 1-D 
hydro simulation such as that used bv IBlondin fc Ellisonl 
l)200l[) . The simu lation is more general tha n the analytic 
approach used bv IDecourchelle et al 1 ll200(t) since it is not 
restricted to self-similar evolution and allows for a con- 
tinuous change in the acceleration efficiency as the SNR 
evolves. We show, however, that when acceleration effi- 
ciency is nearly constant during the self-similar phase, the 
two models closely correspond, providing an important 
check on the validity of both models. While we do not 
calculate X-ray thermal spectra in this paper, we show, 
with various examples, how the efficient production of CR 



protons influences SNR evolution and discuss the implica- 
tions this has on the interpretation of X-ray observations. 



2. CR-Hydro Model 

Our CR-Hydro model couples a spherically symmetric hy- 
drodynamic simulation with a calculation of nonlinear dif- 
fusive shock acceleration. The particle acceleration calcu- 
lation determines how energy is divided between the ther- 
mal gas and relativistic particles, and provides the par- 
ticle distribution over all energies behind the shock, as 
well as the effective ratio of specific heats, 7 e ff, defined 
below. In future work, we will use the electron and ion 
distribution functions to calculate the broad-band contin- 
uum photon emissi on from radio to TeV energies (e.g., 
lEllison et all l200l|) . and use the self-consistent thermal 
prop erties in a non-equilibrium calculation of X-ray lines 
(e.e.. IDecourchelle et all 1200(1) . For now, we restrict our- 
selves to calculating CR proton spectra integrated over 
the SNR lifetime. 

2.1. Hydrodynamic Simulation 

We use a standard hydrodynamic simulation in one dimen- 
sion to model the effects of a s upernova explosion in the 
interstellar medium (ISM) (see IBlondin fc Ellisonl l200ll 
and references therein). We are free to choose arbitrary 
ejecta and ISM mass density profiles, but to facilitate our 
comparisons discussed below, we adopt the par ameters de- 
termined for SN1006 bv lBerezhko et al] l)2002|) . i.e., we as- 
sume a constant density, constant temperature ISM, and 
take the initial ejecta density profile to be p oc r~ n , with 
n = 7 and a co nstant density plateau at small radii (i.e., 
IChevalierLll98l . 2 While the modeling of particular, young 
SNRs depends critically on the ejecta and ISM densities 
(e.s.. IDecourchelle et all 12000). the general character of 
the results we present here are insensitive to these de- 
tails. Specifically, we begin at some time t sta it < 0.1 i c h 
with undisturbed ejecta and ISM separated by a contact 
discontinuity. Here, i c h = R c h/V c h is the characteristic 
age with R ch = [SMcj/iAirpo)} 1 / 3 , V ch = ^2E sn /M cj , and 
po = (1 + 4/He)w p n p o, where E Bn is the explosion energy, 
Mej is the ejecta mass, n p o is the ISM proton number den- 
sity, /ho is the helium to proton number ratio, and m p is 
the proton mass. 3 

The magnetic field, B, is ignored in our hydrodynamic 
model (we implicitly assume that B 2 /8w is small com- 
pared to the thermal pressure), but it is an important 
parameter for the particle acceleration discussed next. 



2 At the start of the simulation, the ejecta temperature is 
assumed to be low enough to be insignificant. 

3 Throughout this paper the subscript (2) implies values 
upstream (downstream) from the shock. 
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2.2. Nonlinear Diffusive Shock Acceleration 

The full details of th e nonlinear acceleration m odel 
used here ar e give n in iBerezhko fc Ellisonl l|l999|i and 
Ellis on et alJ l)2000|) . This is an approximate, algebraic 
model of diffusive shock acceleration containing the es- 
sential physics of NL acceleration, but which parameter- 
izes important properties of the process such as the injec- 
tion efficiency and the maximum energy particles achieve. 
While more com plete models of nonlinear shock acceler- 
ation exist (e.g.. I Jones fc Ellisonl Il99ll IBerezhko et all 
Il99fil | Malk ov fc Dru rv. 2001), our algebraic approxima- 
tion is easier to include in global models of SNRs. It is 
computationally fast making it far less time consuming to 
do parameter searches and to compare model results with 
observations. In Sect. 12.51 we show by direct comparison 
that it gives s imilar results to the m ore physically com- 
plete model of IBerezhko et alJ l)2002l) . 

Briefly, the nonlinear effects in diffusive shock accel- 
eration are: (i) the self-generation of magnetic turbulence 
by counter-streaming energetic particles. Back streaming 
particles produce turbulence in the magnetic field which 
leads to stronger scattering of the particles and hence to 
more acceleration, quickly leading to saturated turbulence 
levels near SB/B ~ 1 in strong shocks; (ii) the modifica- 
tion (i.e., smoothing) of the shock precursor by the back- 
pressure of energetic particles. The precursor influences 
the subshock compression, r su b, the injection and acceler- 
ation efficiencies, and the shape of the accelerated spec- 
trum. Since particle diffusion lengths are generally increas- 
ing functions of momentum (e.g., IBlandford fc Eichleri 
Il987t Idacalone et all Il993h . high momentum particles 
sample a broader portion of the flow velocity profile, and 
hence experience larger effective total compression ratios, 
r to t, than low momentum particles. Consequently, higher 
momentum particles have a flatter power-law index than 
those at lower momenta and can dominate the pressure 
in a NL fashion. The resultant superthermal distribution 
has a characteristic concave upward curvature until the 
spect rum turns over at the highest energies from losses 
(e.g : lEllisori fc Eichleri frQR4l iBlasil l2002UMalkov et all 
2002J); and (iii) the increase in rtot from relativistic par- 
ticle pressure and particle escape. As relativistic parti- 
cles are produced and contribute significantly to the total 
pressure, their softer equation of state makes the shocked 
plasma more compressible (7 — > 4/3). Even more impor- 
tant, as the highest energy particles escape from strong 
shocks they drain away energy flux which must be com- 
pensated for by ramping up the overall compression ratio 
to conserve the fluxes. Just as in radiative shocks, this 
is equivale nt to 7 -> 1 and r tn t. can become arbitrarily 
large Ce.g..lKazanas fc Ellisonl Il98ri IBerezhko fc Ellisonl 
1999t iMalkovt Il997l) . As the overall compression increases 
(''tot > 4), the subshock compression ratio, r su b, which is 
responsible for heating the gas, must become less than the 
test-particle (TP) value (r su b < 4), causing the tempera- 
ture of the shocked gas to drop below TP values. These 
changes in shock compression occur simultaneously with 



changes in the shape of the accelerated particle spectrum, 
thus linking X-ray heating to cosmic-ray p roduction. Fo r 
reviews on diffusive shock acceleration see | Dru ry Jlj^ ) ; 
Blandford fc Eichleri <Tl987h : IBerezhko fc Krvmskvlll 19881) : 



Jones fc Ellisonl I l|l99ll) : and lMalkov fc Drurvl \i 



The most important parameters associated with non- 
linear shock acceleration are the Mach numbers (i.e., the 
shock speed, uq, pre-shock hydrogen number density, n p o, 
and preshock magnetic field, Bq), the injection efficiency, 
r/inj (i.e., the fraction of total protons which end up with 
superthermal energies), and the max imum proton en- 
ergy p roduced, E max . As described in IBerezhko fc Ellisonl 
( 1999), our model includes Alfven heating in the precursor 
which reduces the efficiency compared to adiabatic heat- 
ing and makes the magnetic field strength an important 
parameter. For given sonic and Alfven Mach numbers (i.e., 
given M s = yj p ul / '(7 P ) and M A = y/inp^/ B ), and 
a given shock size and age, Tfrnj sets the overall acceleration 
efficiency and determines the importance of NL effects. 
With other parameters fixed, Alfven wave heating causes 
the acceleration efficiency to decrease with increasing B . 
In a complete model of diffusive shock acceleration, the 
injection efficiency would be determined from first princi- 
ples. However, no current model of diffusive shock accel- 
eration can do this and injection remains dependent on 
approximations of poorly understood wave-particle inter- 
actions. 4 Here, we investigate the effects of efficient accel- 
eration by varying rji n j . A principle aim of future work is 
to constrain ?7i n j from models using X-ray and broad-band 
observations of particular SNRs. 

In order to compa re our results directly to those of 
IBerezhko et all l)2002ll . we assume as they do that the 
magnetic field is turbulent, adopting the Bohm limit for 
strong particle scattering, and somewhat arbitrarily take 
the field downstream from the shock to be the compressed 
upstream magnetic field i.e., Bi = r to tBo, where B0CB2) 
is the upstream (downstream) magnetic field strength. We 
do not expressly consider shock obliquity, i.e., the angle 
between the local shock normal and _Bq, even though this 



4 While in principle plasma simulations, where particles 
move in response to Newton's and Maxwell's equations 
(particle- in-cell) , can fully describe wave-particle interactions 
and injection, in practice, they have not yet done so because 
of computational limits. The limits have been insurmountable 
for two main reasons: (i) The simulations must be performed 
fully in 3D because 1- or 2-D simulations unphysical l y prev ent 
cross-field diffusion (Uokipii et al.lll993l : ljones et al.lil99gfi . In 
all cases except strictly parallel shocks (where the upstream 
magnetic field is parallel to the shock normal), cross-field dif- 
fusion will be an essential part of the injection and acceleration 
process; and (ii) In order for NL effects to become apparent 
or field amplification to occur on large scales, the simulations 
must be run long enough in a large enough box with enough 
particles for a significant population of superthermal particles 
to be produced. If, in addition, electrons are to be investigated, 
the simulations must use the short electron time-step yet run 
for many proton time-scales increasing computation time con- 
siderably. 
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may be an important factor for un derstanding emis sion 
around the rims of some SNRs (see iRevnoldsl ITflfll for 
a discussion of the effects of shock obliquity in a test- 
particle description of particle acceleration in SNRs). As 
a crude approximation, we could model the asymmetry 
seen in many SNRs, including SN 1006, by combining re- 
sults for different quadrants of the remnant where values 
of the magnetic field and injection parameter were varied. 

In our examples presented here, we take the unshocked 
IS M field to be Bp = 20 to match the value determined 
bv lBerezhko et all l|2002|) for SN 1006. For simplicity, un- 
less explicitly stated we use the same constant value for 
the field in the unshocked supernova ejecta even though, 
in reality, this field is likely to weaken considerably with 
time because of flux conservation [in the discussions asso- 
ciated with Figs. H (Sect. ESI and El (Sect. E3, we show 
some effects of a weak ejecta field]. 

The maximum energy cosmic rays obtain depends, in 
part, on the scattering mean free path, A, which is assumed 
to be, 



A = 7], 



mfp Tg j 



(i) 



where rj m { p J> 1 is taken to be a constant and r g = pj (qB) 
is the gyroradius in SI units. Small values of ?7 m fp imply 
strong scattering and allow higher maximum proton ener- 
gies in a given system. The Bohm limit implies rj m f p ~ 1. 

The phase-space mo mentum distributions for protons, 
f(p), are calculated as in lEllison et alJ l 2000h and consist 
of a thermal component, a three-component power law 
at superthermal energies, and a turnover at the highest 
energies given by 



exp 



a \p 



(2) 



Here, a is a constant and p m ax = £ m ax/c is determined by 
setting the acceleration time equal to the SNR age, i sn r, 
or by setting the diffusion length of the highest energy 
particles equal to some fraction, / s k, of the shock radius , 
i? s k ; whichever gives the lowest p max (see iBaring et all 
Il999|) . With these assumpti ons, jtWy i s pro portional to 
the magnetic field strength. IRevnoldsl l|l998j) has shown 
that when synchrotron emission is summed over space in 
a spherically symmetric SNR model, the resultant photon 
spectrum tends to be broader than that produced by a sin- 
gle electron distribution falling off as exp (— p/p maK ). This 
behavior can be approximated in our model by varying a 
and a ~ 0.5 has been used to fit X-ray and radio spectra in 
parti cular SNRs fe.g..lBerezhko et al.lll999l : lEllison et all 
l200l|) . However. Berezhko et all l|2002|) obtain a good fit 
to SN1006 with a sharp turnover in f(p) and to match 
their results, we take a = 4. 

Fig. H illustrates the essential differences in particle 
spectra between TP shock acceleration (dashed line) and 
NL acceleration with different injection efficiencies (solid 
and dotted lines). Compared to the TP power law, i.e., 
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Fig. 1. Schematic particle distribution functions, f(p), 
versus momentum, p, where p • f(p) is plotted to em- 
phasize the spectral curvature in the nonlinear results 
(solid and dotted curves). The kinks in these spectra 
are manifestations of the piec ewise approximation used 
in lBerezhko~ fc Ellison (1999). The lowest momentum 
kink separates the thermal from superthermal population. 
Above that is the three-component power law, and the 
turnover at the highest momentum. The up and down ar- 
row on the test-particle result (dashed curve) indicates 
that there is no constraint in TP models on the density of 
superthermal particles other than that their energy den- 
sity be insignificant. 



the high-energy portion (p > m p c) of a NL spectrum is 
flatter, the low-energy portion of the superthermal spec- 
trum (above the thermal peak and below m p c) is steeper, 
and the thermal part is at a lower temperature. In the TP 
approximation, the normalization of the power-law por- 
tion of the spectrum (relative to the thermal peak) is ar- 
bitrary as long as it is low enough to contain an insignifi- 
cant fraction of the total energy. If we arbitrarily set this 
fraction at 0.01, we can determine, for a given set of shock 
parameters, the injection efficiency, 77tp, below which the 
thermal gas is decoupled from the power law, i.e., for any 
^inj < ^tp, the temperature and density of the thermal 
gas remains constant. For the example shown in Fig. ^ 
r] TP ~ 2.5 x 10 -6 . 

In contrast, the energetic portions of NL spectra are 
set by the conservation of mass, momentum, and energy 
fluxes and any change in the relativistic population influ- 
ences the thermal population. When acceleration is effi- 
cient, the relation between the shock velocity, uq, and the 
temperature of the shocked gas, T2, is no longer approxi- 
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mated by the commonly used expression for a strong, TP 
shock. That is, 

kT * + A , (4) 
fj,m p UQ 16 ' 

where \x m p is the mean particle mass and k is Boltzmann's 
constant. If the equality in equation 0} is assumed for a 
fully ionized plasma with nn e /n p o = 0.1 (njje is the un- 
shockcd number density of helium) and equal downstream 
temperatures for all species, 



1.3 x 10 7 



u 



10 3 kms 



K . 



(5) 



For young SNRs, this predicts temperatures just behind 
the blast wave well above those capable of explaining ob- 
served X-ray line emission (e.g., u — 3000 km s _1 — > 
T 2 ~ 1.2 x 10 s K) and suggests that the shocked elec- 
tron temperature is considerably less than predicted by 
equation jp ) (e.g.jHughes et al I boOOtlHwang et al.l l2002: 
lLong et alTl2003[) . 

Normally, it is assumed that electrons obtain a much 
lower temperature than the ions in the shock layer 
and take some time to equilibrate downstream thro ugh 
Coulomb collisions (see iDecourchelle fc EllisonL l200li for 
a discussion of electron-ion equilibration in modified 
shocks). However, the heating process in the shock layer 
is certainly dominated by collisionless, wave-particle inter- 
actions which are poorly understood. The simplest possi- 
bility, that all particles upon crossing the shock gain a 
speed differential Au ~ uq — U2, where u 2 is the down- 
stream flow speed as seen from the shock, predicts that 
ions are heated far more than electrons. 5 Nevertheless, 
the possibility exists that electrons equilibrate faster than 
this simple conjecture implies and the fact that efficient 
acceleration produces a considerably lower proton temper- 
ature (Fig. 2]) in strong shocks than the TP prediction of 
equation (jSJ), suggests that it may not be as obvious that 
T e2 /Tp 2 <C 1 immediately behind the shock. 6 

Furthermore, radio synchrotron observations have long 
confirmed that many SNRs accelerate electrons to rel- 
ativistic energies, although it is not certain that these 
electrons are drawn from the shock heated population in 
all cases. What is certain is that if relatively cold up- 
stream electrons enter the downstream region without be- 
ing heated and then slowly equilibrate as they advect away 
downstream, very few if any will be injected into the Fermi 

5 That the situation is clearly more complicated than this 
is evident from the fact that if both species only receive 
Am upon crossing th e shock, energy is not conserved (see 
IJones fc EllisonL 1198^ . 

6 Estimates of T e 2/7p2 based on optical and UV line obser- 
vations have been obtained just behind a few SNR shocks (e.g . , 
iRavmond et al.Lll995l : lLaming et allll99dlGhavamiarJ, H999'). 
These show a wide range, i.e., 0.2 < T C 2/T P 2 < 0.8, with large 
uncertainties and are generally restricted to regions that are 
partially neutral where th e most efficient par ticle acceleration 
may not be occurring (see lDrurv et alll200lL for a more com- 
plete discussion). 
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Fig. 2. Parameters from the hydro simulation versus ra- 
dius at various times during the simulation as labeled. 
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from particle acceleration. We have assumed E sn = 3 x 
10 51 erg, M cj 
n p0 = 0.3 cm 

ture is initially a power-law with index n = 7 beyond the 
plateau at small radii (panel c). 



= 1.4 Mq, and a constant density ISM with 
3 , and /hc = 0.1. The ejecta density struc- 



acceleration process. In the simplest, k inematic models 
of injection in Ferm i acce leration (e.g., iJones & Ellison. 
Il99lt iBaring et all [l999), the efficient acceleration of 
thermal electrons is inconsistent with a lack of electron 
heating in the shock layer. 

2.3. CR-Hydrodynamic Coupling 

When efficient particle acceleration shifts energy from the 
thermal gas into relativistic particles, the fraction of to- 
tal energy density in pressure is lowered and the evolution 
of the SNR is modified. The modified evolution, in turn, 
modifies the shock parameters which determine the ac- 
celeration. We model this coupling by including nonlinear 
shock acceleration in a standard 1-D hydrodynamic sim- 
ulation as follows. 

Briefly: (i) the hydro simulation is initiated at some 
time, tstart < 0.1 i c h, after the supernova explosion with a 
plateau-power-law ejecta distribution and a constant in- 
terstellar medium density, pism! (ii) At each time step, 
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Fig. 3. Shock radius versus SNR age for forward (FS) and 
reverse shocks (RS) obtained with our hydrodynamic sim- 
ulation. The dashed curves show NL results with efficient 
particle acceleration (77^ = 10~ 3 ), while the solid curves 
are test-particle (TP) shocks with no particle acceleration. 
The dotted curve shows the TP forward shock radius ob- 
tained with the approxim ate analytic expression given in 
iTruelove fc McKed l|l999|) . 

the hydro evolves and the radii and Mach numbers of the 
forward and reverse shocks are determined; (iii) With this 
shock information, the CR acceleration is calcul ated using 
the approximate model of lEllison et al ] (l200Ctl with the 
injection parameter, ?7i n j, kept constant during the sim- 
ulation and having the same value at the forward and 
reverse shocks; 7 (iv) The acceleration calculation provides 
the overall compression ratio, ftoti an d the effective ratio 
of specific heats is determined by 



7eff 



Mj(n 



2r to 



M|(r to 



1) 



(G) 



where we note that 7 e ff < (P/p cn ) + 1 since j e g includes 
the effects from "escaping" particles. 8 Here, P (p en ) is the 

7 In an actual SNR, of course, the injection efficiency might 
vary with time, vary over the shock surface, or be different at 
the fo rward and reverse shocks (as in our model of Kepler's 
SNR; iDecourchelle et all l2000t) . For the general results pre- 
sented here, keeping Tjjnj constant avoids unnecessary complex- 
ity. It is important to note, however, that the acceleration ef- 
ficiency, i.e., the fraction of ram kinetic energy that is trans- 
ferred to superthermal particles, depends on the shock param- 
eters (Mach numbers, shock age, etc.) as well as on rji n j and 
does vary with age and between the forward and reverse shocks 
even if r;i n j is constant. 

8 Performing the acceleration calculation at each time step is 
computationally expensive. If the full spectrum is not needed 



pressure (energy density) of the shocked gas; (v) The cou- 
pling between the CR acceleration and the hydrodynam- 
ics is accomplished by replacing the ratio of specific heats 
used in the hydrodynamic equations with 7 e g in the shells 
spanning the shock transition. We use a Lagrangian mode 
in the hydro simulation so the simulation grid moves with 
the mass. Once a mass element has been shocked and a 
7eff has been assigned to it using equation J^J, that mass 
element retains that "f e g for the rest of the simulation or 
until it recncounters a shock. Even though energy is leav- 
ing the system via the highest energy cosmic rays, we do 
not explicitly remove energy from the hydro. Lowering 7 e ff 
causes the gas to be more compressible and mimics the es- 
caping energy in one step; (vi) As the SNR evolves, the 
proton distribution functions in momentum, f(j>), asso- 
ciated with each shocked shell undergo adiabatic losses; 
and finally, (vii) At the end of the simulation, the f(p)'s 
from all of the shells are summed to determine the total 
contribution to the Galactic CR source population. 

To start, we show in Fig. [21 a TP case where 7 e g = 
5/3 always (top panel). For this example and all others 
presented in this paper, we have taken E sn = 3 x 10 51 
erg, M e j = 1.4 Mq, Bq = 20 pG, and a constant density 
ISM with n p o = 0.3 cm~ 3 , and /h c = 0.1. These values 
give a characteristic age, t c h = Rch/V c h — 210 yr. The 
initial ejecta density distribution is a power law {p cx r~ 7 ) 
with a flat plateau at small radii. The simulation spans 
6000 yr and various quantities are shown, as a function of 
radius, at different times during the SNR evolution. The 
temperature shown in panel (b) is the average temperature 
defined as 



P U TO„ 



(7) 



where P is the pressure, p is the density, p = (1 + 
4/Hc)/(2 + 3/He)j /He = 0.1 is the ratio of helium atoms to 
protons, and / t h is the fraction of pressure in thermal par- 
ticles. In this plot and all others, we assume the plasma to 
be fully ionized. When acceleration occurs, superthermal 
particles contribute to the total hydrodynamic pressure, 
P. We define the temperature from that fraction of total 
pressure in the "thermal" particles, fth- This fraction is 
determined from f(p) (see Fig. ^ wriere the division be- 
tween thermal and superthermal particles is indicated). 

For the particular parameters assumed for the model 
shown in Fig. |3 the SNR is still in the self-similar phase 
after 100 yr as the reverse shock moves through the power- 
law portion of the ejecta. By 300 years, the reverse shock 
is moving through the plateau region of the ejecta and 
the forward shock has reached Rfs — 4 pc. The density of 
the unshocked ejecta is now well below the ambient ISM 
density Pism = l-4rn p n p o at this stage. By 1000 years, 
Rfs ~ 8 pc, the RS is now moving back toward the origin, 
and the energy of the explosion is shifting from being in 



and only efficient acceleration is considered, a faster approxi- 



mation can use not = 1.3M 3 / 8 if Ma < Ml or r to t = 1.5Aff 



otherwise (see lBerezhko fc Ellisonl 1199' 
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the kinetic energy of the ejecta to thermal pressure. By 
6000 years, the RS has collapsed back through the ejecta 
to the origin and the SNR is well into the Sedov phase 
with a FS speed of about (4/3)750 km s^ 1 . 9 

In Fig. we illustrate how particle acceleration influ- 
ences the SNR shock evolution by comparing TP results 
(Vmj ^ 10~ 6 ) against a simulation with i] in j — 10~ 3 . The 
figure shows the forward and reverse shock radii versus 
remnant age, t snr , for the same parameters as used in 
Fig. El 10 The TP results, i.e., those with j cS = 5/3 every- 
where, are shown with solid curves and the NL results are 
shown with dashed curves. The dotted curve is the TP an- 
alytic approximation given bv iTruelove fc McKeei l|l999|) 
for the FS and this matches the hydro results well over 
the entire time period. When NL effects are important, 
the region between the forward and reverse shocks be- 
comes narrower and d enser, a result fully consis tent with 
the analytic results of iDecourchelle et al. I ll2000h and the 
multi- dimensional hydro simulations of lBlondin fc Ellison 
lj200lh . 

In Fig. 21 we compare TP profiles to NL ones at 
i snr = 200 and 1000 yr. In all panels, the heavy-weight 
solid curves are the TP results (rfinj ^ 10~ 6 ), the dotted 
curves have ?7inj = 1-2 x 10~ 4 , and the dashed curves have 
?7i n j = 10 -3 . Panels (a) and (e) show the reduction in 7 e g 
that occurs as shock acceleration shifts energy into rela- 
tivistic particles. This produces compression ratios greater 
than 4 and results in denser, narrower regions between the 
forward and reverse shocks as shown in the density pro- 
files (6) and (/). The separation between the forward and 
reverse shocks varies considerably with i snr and rji n j. At 
t snl = 200 yr, this separation is ~ 0.25 pc for = 10~ 3 
and ~ 0.7 pc for the TP case, while at t snl — 1000 yr it 
is ~ 4.2 pc for 77^ = 10~ 3 and ~ 6 pc for the TP case. 
The shift in energy out of the thermal gas also results in a 
significant decrease in the temperature of the shocked gas 
in the NL shown in panels (c) and (g). 

The flow speed is shown in the bottom panel and the 
forward and reverse shocks are easily identified in all pan- 
els. At t snl — 1000 yr, the flow speeds just behind the for- 
ward shocks of the three examples are about equal, i.e., 
Ud s (TP) ~ 2600 km s" 1 and u ds (?7inj = 1CT 3 ) ^ 2500 km 
s" 1 . However, since the shock speed depends on r to t, the 
TP shock is moving faster, i.e., m (TP) ~ (4/3)2600 ~ 
3500 km s-\ while u (mn S = 10~ 3 ) ~ (7/6)2500 ~ 2900 
km s . 

The intermediate NL case with j^j = 1.2 x 10~ 4 
(dotted curves) produces temperatures and compression 

9 Note that the speed shown in the bottom panel of Fig. |5] 
and panel (d) of Fig. 2]is the flow speed in the frame of the ex- 
plosion. The shock speed is given by Mo = (notiids— M up )/(rtot — 
1), where u up (itd s ) is the flow speed just upstream (down- 
stream) from the shock and r to t is the shock compression ratio. 
For the ISM conditions we assume here, M up always equals zero 
for the forward shock. At t Bnr = 6000 yr in Fig. |5]at the FS, 
uds — 750 km s~\ r to t — 4, and Mo m 1000 km s _1 . 
10 We do not follow particle acceleration at the reverse shock 
after it first collapses to zero radius. 



ratios between the TP and ^; n j = 10~ 3 cases, as ex- 
pected. However, as seen in the variations of j c g (pan- 
els a and e of Fig. 0}, the ?/ in j = 1.2 x 10~ 4 case acts 
quite differ ently at the beginning of the simulation. As 
described in lBerezhko Ez Ellison! l|l999f) . unmodified shocks 
with r to t — 4 can occur for high Mach numbers if n- m j is 
small enough. In this case, the pressure in relativistic par- 
ticles, P rc i, is small compared to PqUq and these particles 
do not slow the incoming gas enough to produce the non- 
linear shock modification. For a given ?yi n j, as the SNR 
ages, Mo decreases, P r ei/(po«o) increases, and the shock 
acceleration becomes more efficient and more nonlinear. 
Thus, for relatively inefficient injection, 7 e ff drops and r to t 
increases in the early stages of evolution, exactly the op- 
posite behavior as with efficient injection. 11 

After the initial transition from a strong, unmodified 
shock to a strong, modified one, the forward shock will 
enter the long stage where it slowly weakens, the accel- 
eration becomes less efficient, and 7 c ff increases toward 
5/3. However, as shown by the dotted curve in panel (e) 
in Fig. 01 the reverse shock weakens much more quickly 
and at t snl = 1000 yr has 7 G fr — 5/3. Since we assume 
that Bq is a constant everywhere, the RS weakens as the 
ejecta density drops and the magnetic pressure becomes 
dominant over the gas pressure. In an actual SNR, mag- 
netic flux conservation in the expanding, unshocked ejecta 
would cause B to weaken as the ejecta density drops and 
the RS would remain stronger than with a constant Bq. 
For a given unshocked density, a weak Bo results in effi- 
cient shock acceleration because the Alfven Mach number 
increases with decreasing -Bo, and because the amount of 
energy transferred out of accelerated particles into waves 
and then into heating of the precursor gets less. Of course, 
if the magnetic field gets too weak, the gyroradii of ac- 
celerated particles will become comparable to the shock 
radius, or the acceleration time will become comparable 
to t snr , before high particle energies or high efficiencies are 
obtained. To illustrate the effects of a low ejecta magnetic 
field, we show in panels (a) and (e) of Fig.^Heff produced 
by our rj- m ^ = 1.2 x 10 example with the unshocked ejecta 
field set at 1 pG (light-weight solid curves). As seen clearly 
in panel (e), j c s remains lower at the end of the simula- 
tion than in the (-Bo)rs = 20 pG case, indicating more 
efficient acceleration. 

In Fig. we plot the compression ratio, r to t, versus 
t snl for our TP and two NL examples. There are two sets 
of curves because there are two ways of determining r to t ■ 
The solid curves are values obtained from the particle ac- 
celeration calculation which are then used in equation (jSJ 
to calculate 7 c g. The dashed curves are determined di- 



11 As discussed in iBerezhko fc Ellison! (|T999), it is by no 
means certain that actual SNR shocks will have injection rates 
low enough for high Mach number, unmodified shocks to occur. 
Furthermore, there may be large differences in injection rates 
for parallel and oblique regions of the same SNR blast wave, 
making it possible that some regions are highly modified and 
others unmodified. 
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Fig. 4. Comparison of the parameters for a test-particle simulation (solid curves) with those undergoing particle 
acceleration. The results are shown at t sal — 200 and 1000 yr and the supernova parameters are the same as those 
in Fig. The light-weight solid curves in panels (a) and (e) show the effect on j c g of reducing the unshocked ejecta 
magnetic field from (-Bo)rs = 20 /iG to 1 /iG. Apart from (-Bo)rSj all other parameters are the same as in the 
77 in j = 1.2 x 10 -4 examples shown with dotted curves. In panel (e), the shocked ISM portion of the (£>o)rs = 1 /iG 
curve is essentially identical to the dotted curve. 



rectly from the hydro density ratio taken just downstream 
and upstream from the shock, /Cdw/pup- In principle, these 
should be identical but in practice differences occur, indi- 
cating the inherent errors in our 7 e ff technique. The largest 
average difference shown here with ?7i n j = 10~ 3 is < 10%. 
The ?7i n j = 1.2 x 10~ 4 case clearly shows the rapid tran- 
sition from the strong, unmodified shock with r to t ~ 4 
at early times to the modified shock with r tot > 4 at 
intermediate times. The striking difference between the 
??inj = 10~ 3 and ?7i n j = 1.2 x 10 -4 cases will leave an im- 
print on the post-shock gas that may be an important 
diagnostic for determining if efficient Fermi acceleration 
takes place in young SNRs. 



In Fig. we show the postshock proton temperatures 
for the same three examples. Comparing the heavy-weight 
solid and dashed curves shows that the RS temperature in 
the ?7i n j = 10~ 3 efficient shock is less than ~ 1/10 that of 
the test-particle shock over most of the time period shown. 
This is a difference large enough to have profound conse- 
quences for the interpretation of X-ray thermal emission 
in young SNRs. Furthermore, the temperature in the in- 
termediate case (r^inj = 1.2 x 10 -4 ) shows a stronger time 
variation than either the TP or r]i n j = 10 -3 case. This is 
reflected as a somewhat steeper spatial gradient in tem- 
perature as seen in panel (g) of Fig. 0] again offering a 



Donald C. Ellison, Anne Decourchelle and 




. Solid: Accel, cal 
. Dashed: Hydro 

Q I l i I i i I i i l l l 

10 100 1000 

Fig. 5. Forward shock compression ratios vs. supernova 
age for various injection rates, ?7inj. In all cases, the solid 
curves are obtained from the NL particle acceleration cal- 
culation and the dashed curves are obtained from the den- 
sity ratio, Pdw/Pupj just downstream and upstream from 
the shock. 
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Fig. 6. Temperatures immediately behind the forward 
and reverse shocks as a function of i snr . The solid curves 
are the TP case, the dotted curves are for ?7i n j = 1.2 x 10~ 4 , 
and the dashed curves are for ?7i n j = 10~ 3 . The reverse 
shock temperatures are plotted with heavy-weight curves. 



possible diagnostic to test for efficient particle accelera- 
tion. 

2.3.1. Comparison with Analytic Model of 
Decourchelle et al. 

In our previous work l|Decourchelle et all l200o|) we cou- 
pled NL acceleration with an analyti c, self-similar des crip- 
tion of the SNR hydrodynamics (i.e. . IChevalierL fl983|) and 
a non-equilibrium ionization calculation of thermal X-ray 
emission. That model used the same algebraic acceleration 
calculation we use here with the same injection parame- 
ter, ?7i n j. The assumptions of the two models differ mainly 
in the way the hydrodynamics are treated. 

While the hydro simulation has a one-fluid approach 
with an effective ratio of specific heats derived from the 
compression ratio (equation EJ, the analytical, self-similar 
solutions are based on a two- fluid approach: a nonrelativis- 
tic one with 7 = 5/3 and a relativistic one with 7 = 4/3. 
The NL particle acceleration calculation provides the com- 
pression ratio and the fraction of total pressure in relativis- 
tic gas at both shocks which are used to determine the 
boundary conditions of the self-similar solution. Imposing 
the compression ratios, which can be greater than 7, mim- 
ics the effects of particle loss. The limitations of the self- 
similar calculation are that the RS must be propagating 
in the power-law portion of the ejecta density profile, and 
that it assumes that the forward and reverse shock com- 



pression ratios, as well as the fractions of total pressure 
in relativistic gas, remain constant throughout the evolu- 
tion. A limitation of our hydrodynamic simulation is that 
it assigns at the shocks to an element of gas an effective 
gamma, which is kept constant during its post-shock evo- 
lution. 

In the top three panels of Fig.Qwe compare the density 
structures for our two models at t snr = 100 yr when the RS 
is still propagating in the power-law portion of the ejecta 
density profile. We use the same supernova and ISM pa- 
rameters as in our previous runs. For the TP (panel a) and 
?7i n j = 10~ 3 (panel b) cases, the correspondence between 
the two models is extremely close with small differences 
resulting mainly from the finite grid used in the hydro 
simulation. The reason for this comes from the fact that 
self-similarity is satisfied in the TP case where r to t — 4 al- 
ways. It is also approximately satisfied in the ?7i n j = 10 -3 
case since r to t, and therefore 7 c fr, remains fairly constant 
during the 100 yr evolution, as indicated by the dashed 
curve in panel (a) of Fig. 

While it was not obvious that the 7 e fi procedure used 
in the hydro model would be a good approximation, the 
excellent agreement we see with the analytic model, which 
only assumes self-similarity, justifies this approach, at 
least at these early times. 

Much larger differences in the two models are seen in 
panels (c) and (d) where ?7i n j = 2 x 10~ 4 . In this case, 
rtot and therefore j e s vary strongly during 100 yr (similar 
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Radius [pc] 

Fig. 7. Density (top three panels) versus radius for two 
NL models and a TP model, as indicated, at i snr = 100 
yr. The bottom panel shows temperature versus radius for 
?7i n j = 2 x 10~ 4 . In all panels, the solid c urves are calculated 
using the analytic, self-similar model of lDecourchelle et alJ 
(2000) and the dotted curves are from the CR-Hydro 
model presented here. In the CR-Hydro model, some nu- 
merical noise is evident in the shocked ISM profiles for 
r, inj = 2 x 10- 4 . 

to the dotted curve in panel a of Fig. , and self-similar 
conditions no longer apply. Since the analytic model takes 
r to t at the forward and reverse shocks at the final time and 
applies these values for the entire evolution, it is effectively 
using values, at least for the FS in the 77j n j = 2 x 10~ 4 case, 
which are considerably greater than the average values 
used in the hydro calculation. Using a larger rtot means 
there is less pressure for a given energy density pushing the 
shock in the self-similar model than in the hydro model. 
Therefore, the FS in the hydro calculation travels faster 
and extends further than in the analytic calculation. The 
variation in 7 e g is less in the shocked ejecta (again similar 
to the dotted curve in panel a of Fig. so the difference 
in the two models is less there as well. 

We consider the close correspondence between these 
two models in the TP and 7yi n j = 10~ 3 cases to be a clear 
indication that both adequately model the SNR hydro- 



dynamics during self-similar conditions. The validation 
of the 7 c ff approach by our two-fluid description of the 
post-shock flow in the framework of self-similar, cosmic- 
ray modified Chevalier solutions means the more general 
hydro simulation can be used with some confidence when 
self-similar conditions do not apply. This will also allow us 
to use the analytic model to test the hydro model when X- 
ray thermal and broad-band photon emission is included. 

2.4. Energetic Proton Spectra 

The proton distribution functions that are calculated at 
the forward and reverse shocks during the SNR evolution 
can be summed to determine the total contribution to the 
cosmic-ray spectrum. The amount of material swept up by 
each shock is used to weight each spectrum and convert 
f(p) to N(p,At) where 4tt p 2 N(p, At)dp is the total 
number of particles overtaken by the shock in the time in- 
terval, At. The iV(p, At)'s are then added together, with 
each one adjusted for adiabatic losses during the time from 
when the spe ctrum was calcu lated until the end of the sim- 
ulation (e.g.. iRevnoldsl fl99^1 . to produce N(p). 12 Typical 
proton spectra calculated with Bq = 20 /iG at t snI — 500 
yr (WU =i 2.4 for E sn = 3 x 10 51 erg, M ej = 1.4 M , 
n p0 = 0.3 cm~ 3 , and /h = 0.1) are shown in Fig. |H1 
The NL examples produce cosmic rays at the expense of 
heating and show lower temperatures, as indicated by the 
positions of the thermal peaks. The dashed curve is the re- 
sult with the most efficient injection (rft n j = 10 -3 ) and this 
spectrum is the flattest at relativistic energies (before the 
turnover at p max ). The spectrum with ?7i n j = 1.2 x 10 -4 
is the most distorted from a Maxwellian at thermal en- 
ergies due to the summing of spectra which evolve from 
unmodified (rtot — 4) to strongly modified (rtot > 4), as 
discussed above. 

To illustrate the effects of a weakening unshocked 
ejecta magnetic field, we show in the bottom panel of 
Fig. El (thin solid curve) N(p) for the case where the un- 
shocked ejecta field is set at (Bo)rs = 1 A*G, a factor 
of 20 lower than in the other cases. The unshocked ISM 
field is kept at 20 /j,G and ?7i n j = 1.2 x 10~ 4 for both 
the forward and reverse shocks. The flatness of the su- 
pcrthermal part of the spectrum compared to the dotted 
curve indicates that the RS is stronger (i.e., r to t is larger) 
than when (-Bo)rS = 20 fiG (see Fig.0J. However, p max is 
well below that in the (-Bo)rs = 20 /iG cases because the 
weak field produces large gyroradii and high momentum 
particles cannot be accelerated. The FS spectrum is only 
slightly changed by changing (£?o)rs and is not shown. 

In Fig. we show spectra at 1000, 6000, and 4 x 10 4 
yr for the same ISM and SN parameters with 77^ = 10~ 3 . 
As the FS expands and weakens, the acceleration efficiency 
drops and the relativistic portion of the spectrum steep- 
ens. As expected, particles accelerated at the RS decrease 

12 We do not include radiation (i.e., synchrotron) losses here 
since they are always negligible for ions. For electrons, of 
course, these losses must be considered. 
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Fig. 8. Proton number distributions, N(p) (multiplied by 
[p/(m p c)] to p /Mq), for hydro simulations with no accel- 
eration (heavy solid curves), with ?7i n j = 1.2 x 10~ 4 (dot- 
ted and thin solid curves), and with r]i n j = 10 -3 (dashed 
curves) summed over 500 years. The spectra in the top 
panel are from the shocked ISM (forward shock), while 
those in the bottom panel are from the shocked ejecta 
(reverse shock). The N(p) , s are absolutely normalized to 
the total number of protons overtaken by the shock during 
their lifetime. All spectra have been adjusted for adiabatic 
expansion. The thin solid curve in the bottom panel shows 
the RS spectrum for the case where the unshocked ejecta 
magnetic field is set at (-Bq)rs = 1 £*G- 



in importance in later stages of the SNR. As we show 
in Fig. II II below, even moderately efficient injection (i.e., 
Vinj 2i 10~ 4 ) results in approximately 50% of the explo- 
sion energy being put into relativistic particles over the 
lifetime of the SNR. Another important property shown 
in the top panel of Fig. is that, even with rji n j = 10 -3 
which produces strongly curved spectra at early times, the 
total particle distribution function after 4 x 10 4 yr is ap- 
proximately N(p) oc p~ 2 at relativistic energies. Thus the 
total cosmic-ray production is fairly consistent with the 
observed CR flux after energy dependent escape from the 
galaxy is included. 
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Fig. 9. Particle number distributions, N(p) (multiplied 
by [p/(m p c)] 4 TO p /M ) from the shocked ISM (FS) and 
shocked ejecta (RS) at various times during the SNR evo- 
lution. Efficient particle injection is assumed (jjmj = 10~ 3 ) 
throughout the SNR evolution, but the fraction of energy 
in relativistic particles varies as the SNR ages. The ac- 
celeration efficiency decreases as the SNR ages and the 
shocks weaken. The RS contributes little to the total CR 
population at t snI = 4 x 10 4 yr. All spectra have been 
adjusted for adiabatic expansion. 



2.5. Comparison with Kinetic Models 

One of the most complete models of nonlinear CR pro- 
duction in SNRs is that of Berezhko and co-wo rkers (e.g., 
iBerezhko etafl Il99fil IBerezhko fc Volkl Il997l) where the 
time-dependent CR transport equations are solved to- 
gether with the gas-dynamic equations in spherical sym- 
metry. This model gives the radial distribution of gas and 
accelerated CR spectrum at any phase of the SNR evo- 
lution and has been used with good suc cess to model a 
numb er of young SNR s including SN1006 feerezhko et"afl 
l2002l). Cass i opeia A l|Berezhko et all l2003j) . and Tvcho 
(Vo lk et alll2002|) . In all of the above cases, the authors 
conclude that satisfactory fits to the radio and X-ray syn- 
chrotron emission, along with 7-rays where observed, are 
best obtained when the forward shocks are strongly modi- 
fied by the efficient acceleration of protons in a high mag- 
netic field. While similar resul t s have been presented be- 
fore (e.g., iRevnolds fc EllisonL Il99l iRlliaon et all EM 
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Fig. 10. The dashed curve is the i nferred FS proton dis- 
tribution of SN1006 as obtained bv lBerezhko et alJ l|2002|) 
at t sm — 1000 yr. Using their input parameters for E sn , 
M e j, etc., we have matched this result (solid curve) with 
?7inj = 1.2 x 10~ 4 , a = 4, and / s k = 0.1 (see equationEland 
the discussion following it). The dotted curve is the cor- 
responding proton distribution from shocked ejecta (RS) 
from our model. 



2001), Berezhko's model probably contains the most com- 
plete physical description of the NL acceleration process 
at the FS, and the detailed modeling of particular SNRs 
clearly lends credibility to the suggestion that shocks in 
young SNRs can be strongly nonlinear with compression 
ratios greater than 4. 

Our CR-hydro model is complementary to that devel- 
oped by Berezhko and co-workers in that we use a more 
approximate model of the NL acceleration, but have a 
more complete model of the hydrodynamics including the 
reverse shock. In particular, we do not follow the particles 
through the acceleration process where previously accel- 
erated particles can continue to be re-accelerated at later 
times. Instead, we assume that all accelerated particles 
originate as thermal particles in the upstream flow, are ac- 
celerated as the shock overtakes them, and then remain in 
the shell of material where they were accelerated, suffering 
adiabatic losses as the remnant ages. The re-acceleration 
of high energy particles with long diffusion lengths can be 
important as it adds to the max imum energy particles ob- 
tain (e.g. iBerezhko et all Il996|) . but this effect becomes 
less important as time goes on and new unshocked parti- 
cles are overtaken and accelerated. 

In Fig. 1101 we c ompare our results with those of 
IBerezhko et al.l ( 2002) for the proton distribution at t sm = 
1000 yr they predict for SN1006. We have adjusted the 



normalization of the FS to obtain a good fit, but oth- 
erwise have used the environmental parameters given in 
Berezhko et al. While both our models require an arbitrary 
injection parameter defined as the fraction of thermal par- 
ticles overtaken by the shock that become superthermal, 
the actual implementation of this parameter may be dif- 
ferent in the two models. To fit SN 1006, Berezhko et al. 
use r\ — 2 x 10~ 4 ; we have used Tftnj = 1.2 x 10~ 4 to 
match their results. A distinct advantage of Berezhko's 
model is that the maximum CR energy, E m&x , and the 
shape of f(p) near E max are determined self-consistently. 
In our CR-hydro model, we parameterize these with a in 
equation for the shape and with the fraction of shock 
radius, / s k, set equal to the maximum particle diffusion 
length to determine E maxi . To match the results given by 
Berezhko et al . (2002), we use a = 4 and f s \. = 0.1. While 
fine tuning of our parameters could provide a more exact 
match than shown in Fig. 1101 the particular values are 
unimportant compared to the fact that both models show 
a similar level of injection efficiency and shock modifica- 
tion. Considering the fundamental differences in the two 
models, the fit of superthermal particles accelerated by 
the FS shown in Fig.^|is impressive. The concave shape, 
accentuated by plotting p 4 N(p), is very similar, although 
the piecewise nature of our approximation is apparent. 

Equally important, the FS compression ratios are very 
similar, as shown in the top panel of Fig. 1111 where the 
total and subshock compression ratios are compared for 
this SN1006 model (CR-hydro: solid curves; Berezhko et 
al.: dashed curves). As we discussed earlier, unmodified 
strong shocks can exist when ?yi n j is relatively low and this 
is just what is predicted in both our models for i snr 200 
yr in SN 1006. As the FS slows, the acceleration becomes 
more efficient, r to t increases, and r su b starts to drop below 
4 at t sm <; 200 yr, indicating that the shocked proton 
temperature will be less than predicted in the TP case. 
The dotted curves are CR-hydro results with various ryinj 's 
and show (curves b and c) that the transition between 
unmodified and modified shocks can be extremely abrupt. 
For large enough r)i n j (curve d) unmodified solutions do 
not occur and r to t is large from the start. 

We take the excellent correspondence achieved with 
Berezhko's calculation as evidence that our CR-hydro 
model, with its algebraic NL particle acceleration approx- 
imation, is accurate. 

As is clear in Fig. B erezhko et all ((20021) do not 
treat thermal particles explicitly although, in principle, 
they could convert the thermal pressure and density into 
a Maxwellian distribution as we do. They also ignore the 
reverse shock, an approximation that is well justified if 
only continuum emission from relativistic particles is con- 
sidered at times t snr 3> £ c h- In order to predict thermal 
X-ray emission in young SNRs, however, thermal elec- 
trons must be modeled and particle heating and accel- 
eration at the forward and reverse shocks must be self- 
consistently determined. We do not show electron spectra 
here, but our NL acceleration model, with additional pa- 
rameters, can produce these along with the broad-band 
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Fig. 11. Forward shock compression ratios (top panel) 
and fraction of E sn going into relativistic particles, 
E CT /E sn , (bottom panel) for SNRs with various r/inj 's 
versus t snr . The solid curves are for SN1006 parame- 
ters (ri in j = 1.2 x 1(T 4 , B = 20 ^G, E sn = 3 x 10 51 
erg, M ej = 1.4M Q , n p0 = 0.3 cm" 3 , and / Hc = 0.1) 
and show r to t and r su b in the top panel and the total 
E cr /E sn and that portion from the RS in the bottom 
panel. The dashed curves are the corresponding SN1006 
results from lBerezhko et all (|2002^) . The dotted curves are 
from the CR-hydro model with (a) 77^ = 8 x 10~ 5 , (b) 
r, inj = 2.4 x lO" 4 , (c) 7? inj = 4 x 10" 4 , and (d) rfoj = 10~ 3 . 



continuum em ission from radio to TeV 7-rays produced 
by them (e.g., lEllison et all EoOlh . The most important 
advantage our CR-hydro model has, however, is the abil- 
ity to model the RS and the potential for including X- 
ray thermal emission from the shock-heated ejecta calcu- 
lated using a non-equilib rium ionization calculation (e.g., 
InecourchelleeirnEoOri. 

In the bottom panel of Fig. ^JJ we show the fraction 
of explosion energy in relativistic particles, E CT /E an , as 
a function of i snr . As in the top panel, the solid curves 
are from our CR-hydro mode l and t he dashed curve (FS 
only) is from IBerezhko et alJ l|2002ft . The total E cr /E sn 
is in good agreement at early times but diverges some- 
what at later times, a possible indication that adiabatic 
losses are treated differently in the two models. Most im- 



portantly, both models show that approximately 50% of 
the explosion energy is placed into cosmic rays during the 
SNR lifetime. 13 The dotted curves in the bottom panel 
show that E cr /E sn at £ snr = 6000 yr remains in a fairly 
narrow range (~ 40 to 65%) for 77j n j ranging from 8 x 10~ 5 
(a), to 10~ 3 (d). 



3. CONCLUSIONS 

The efficient production of cosmic rays by shocks in SNRs 
lowers the pressure to energy density ratio in the post- 
shock gas causing dramatic changes in the thermal prop- 
erties of the shocked gas and the SNR evolution. We have 
presented a CR-hydro model that combines a 1-D hydro- 
dynamic simulation of a SNR, including the forward and 
reverse shocks, with particle acceleration. We explicitly 
include the effects of particle acceleration on the shock 
heated ejecta, a critical step in determining how X-ray 
thermal emission from the hot ejecta is modified by par- 
ticle acceleration. 

In accord with previou s results (e.g., iDorfiL Il990t 
IBerezhko etafl Il996l 12002^1 . we find that SNRs can eas- 
ily transfer ~ 50% of the explosion energy into relativis- 
tic particles (bottom panel of Fig. 1 1 1 ft . Compared to the 
situation where acceleration is absent or inefficient, this 
transfer results in lower FS speeds, larger overall compres- 
sion ratios, cooler post-shock temperatures, and a smaller 
and denser interaction region between the forward and re- 
verse shocks. At young ages, compression ratios r tot > 6 
(top panel of Fig. ITT)) and proton temperatures less than 
£5 1/10 the TP value are predicted (Fig. EJ) for the par- 
ticular set of SN and ISM parameters used here to model 
SN 1006. For some injection efficiencies, an abrupt tran- 
sition is predicted to occur at early times (i snr < t c b) be- 
tween high Mach number, unmodified shocks with r to t ~ 4 
and strongly modified shocks with r tot ^S> 4 (top panel 
of Fig. Ill|) . causing steep spatial gradients in tempera- 
ture (Fig. We believe the changes in shock speed, den- 
sity, temperature, and spatial extent and profile are large 
enough to (i) provide diagnostics for X-ray observations of 
young SNRs sufficient to place meaningful constraints on 
the acceleration efficiency, and (ii) to importantly modify 
the inferred values of E sn , M e j, and ambient ISM density 
once a full thermal X-ray emission model is combined with 
the CR-hydro model. 

Despite the strong NL effects expected for young 
SNRs, the CR spectrum integrated over the age of a rem- 
nant (Fig. |5J| should have a spectrum not much flatter 
than N(E) oc E~ 2 at relativistic energies. 

Our one-dimensional CR-hydro model uses a compu- 
tationally fast, approximate calculation of particle accel- 



13 IBerezhko et alJ J2002I) make the important point that injec- 
tion may vary over the surface of the SNR and be significantly 
less where the magnetic field is highly oblique. They estimate 
that to supply the galactic cosmic-ray population the overall 
efficiency need only be ~ 20% of the maximum values shown 
in Pig. El 
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eration which is coupled to the hydrodynamics by mod- 
ifying the effective ratio of specific heats. We have veri- 
fied the accuracy of this approach for parameters applica- 
ble to SN1006 by direct comparison with the more phys- 
ically complete model of acceleration of iBerezhko et alJ 
l|2002() (Figs. Un| and UTJ, where the time-dependent 
CR transport equations are solved self-consistently with 
the gas-dynamic equations. The main advantage of our 
model lies in the fact that we can model acceleration 
at the forward and reverse shocks during all stages of 
the SNR evolution. In a preliminary work we modeled 
thermal X-ray emission from Kepler's SNR with a two- 
fluid, self-similar description of the SNR hydrodynam- 
ics co upled to the same calcul ation of particle acceler- 
ation l|Decourchelle et all Eooof) . We have demonstrated 
(Fig. EJ that the two models give essentially identical re- 
sults when self-similar conditions apply. Our next step will 
be to include this non-equilibrium calculation of thermal 
X-ray emission, plus broad-band continuum emission from 
bremss trahlung, synchrotr on, inverse-Compton, and pion- 
decay l|Baring et al ., 1999), in the CR-hydro model. 
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